Solving the Fokker-Planck kinetic equation on a lattice 
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O^l , We propose a discrete lattice version of the Fokker-Planck kinetic equation along lines similar to 

the Lattice-Boltzmann scheme. Our work extends an earlier one-dimensional formulation to arbi- 
trary spatial dimension D. A generalized Hermite-Gauss procedure is used to construct a discretized 
kinetic equation and a Chapman- Enskog expansion is applied to adapt the scheme so as to correctly 
. reproduce the macroscopic continuum equations. The stability of the algorithm with respect to the 

finite time-step At is characterized by the eigenvalues of the collision matrix. A heuristic second- 
order algorithm in At is applied to investigate the time evolution of the distribution function of 
simple model systems, and compared to known analytical solutions. Preliminary investigations of 
sedimenting Brownian particles subjected to an orthogonal centrifugal force illustrate the numerical 
efficiency of the Lattice-Fokker-Planck algorithm to simulate non-trivial situations. Interactions 
between Brownian particles may be accounted for by adding a standard BGK collision operator to 
the discretized Fokker-Planck kernel. 
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O ! I. INTRODUCTION 



Kinetic equations are well established mathematical models for investigating the out of equilibrium behaviour of 
fluids, and their relaxation towards thermodynamic equilibrium, at a molecular or coarse-grained, mesoscopic level 0. 
They govern the time evolution of the single particle distribution function /(x, v;t) in the 2 x D-dimensional space 
of position x and velocity v. This evolution is expressed in terms of "free flow" , under the action of an external 
[ or self-consistent force field, and of the action of a "collision" operator C[f], which accounts for the interactions 
t-H ' between particles, or their coupling to a continuous medium. The exact form of the operator C involves a hierarchy 
of equations for the higher order distribution functions (the BBGKY hierarchy @), so that a closed equation for / 
cannot be obtained. Depending on the physical problem at hand, approximate closures have been devised which lead 
to various standard kinetic equations. 

Thus, if particle interactions are only considered at a mean-field level, through a self-consistent force field, (7 = 
and the standard Vlasov equation of plasma physics results [j| . In dilute gases of molecules interacting through short- 
range forces, one may make the assumption of strictly binary, uncorrelated collisions, which leads to the non-linear 
Boltzmann collision operator involving the molecular scattering cross-section Jl| . The Boltzmann equation has been 
widely used for a systematic investigation of transport phenomena in gases while its generalization by Enskog, 
O ■ which accounts for static correlations, allows such calculations to be extended to dense fluids 

Following the idea that the molecular details included in the Boltzmann and Enskog collision operators are not 
likely to have a strong influence on the experimentally measured macroscopic properties of fluids, Bhatnagar, Gross 
and Krook (BGK) proposed a highly simplified, phenomenological version of the collision operator describing the 
relaxation of the distribution function towards local Maxwellian equilibrium on a single time scale r. The BGK 
operator Q still conserves mass, momentum and kinetic energy, and by properly adjusting the relaxation time r, it 
goes beyond the strictly binary collision assumption of the Boltzmann equation and hence allows its phenomenological 
extension to dense fluids 6]. The combination of the BGK kernel with discretized lattice versions of kinetic equations, 
globally referred to as the Lattice Boltzmann (LB) method nas proved to be a powerful tool for the study of 

laminar or turbulent fluid flow and transport. The assumption that fluid particles can be restricted to have only a 
small, fixed number of velocities v reduces the computationtal problem considerably compared to corresponding finite 
element schemes. The numerical parameters of the Lattice-BGK model can be adjusted to reproduce the correct 
Navier-Stokes behaviour in the small velocity (or Mach number) limit. Over the last decade the LB method has 
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increased in popularity as successful applications have been repeatedly reported in numerical simulations of large 
scale hydrodynamic flows 0, complex fluids under shear and in porous media |10| , ion transport in nanochannels , 
and colloidal suspensions 

The latter systems generally involve a considerable separation in size and time scales as epitomized by the classic 
concept of Brownian motion. This is generally the case of two-component systems involving a molecular-scale solvent 
and larger, heavier solutes. The natural kinetic theory framework to handle such highly asymmetric situations is 
the Fokker-Planck (or Kramers) equation [l4|, which adopts an effective, one-component description of the solute, 
whereby the solvent and the boundaries are modelled implicitly, as sources of friction and random forces. The collision 
operator C[f] may then be constructed as the sum of a Fokker-Planck (FP) operator, which accounts for the coupling 
between the solute and the (continuous) solvent and between the solute and the confining surfaces, and of a BGK 
operator to model solute/solute interactions. The corresponding lattice Fokker-Planck (LFP) equation was recently 
put forward by some of us |15| . and applied to a simple one-dimensional problem of electrical conduction. 

The main objective of the present paper is to extend the LFP formulation to the D-dimensional case, and to 
develop an efficient and stable numerical scheme for its solution. This should provide an operational tool to tackle 
non-equilibrium problems in the field of dispersions and complex fluids involving multiple length and time scales. 

The paper is organized as follows. The lattice discretization of the FP equation is carried out in sec. [HJ Using a 
truncated expansion of the distribution function /(x, v;£) in generalized Hermite polynomials, the collision operator 
C[f] is expressed in terms of the moments of /, which can be computed by appropriate quadratures. This completely 
defines the LFP numerical solution scheme. In sec. Mil we address the stability of such a scheme. Since the evolution of 
the discretized distribution functions can be rewritten as a linear iteration, studying the stability amounts to analysing 
the spectrum of the transformation. A standard Chapman-Enskog expansion of the LFP equation is carried out in 
sec. II VI to ascertain the reproducibility of the continuous macroscopic equations. The practical implementation of a 
second order algorithm in the discrete time step At is proposed in sec. while numerical results are presented in 
sec. I VII Concluding remarks are contained in sec IVIII and mathematical issues are detailed in the appendices. 



II. THE LATTICE FOKKER-PLANCK EQUATION 

Since the implementation of the BGK collision operator in the LB method is amply documented in the literature 0, 
0, we restrict the following to the Fokker-Planck operator. The standard FP kinetic equation in D dimensions 
reads Q: 

(d t + v a d a + a E d Va )f = C FP [f] = jd Va (v a + v 2 T d Va )f (1) 

where x a and v a are the cartesian components of the D-dimensional position and velocity vectors, d a and d Va 
the corresponding gradient operators, and a E are the components of the acceleration due to an external force field 
ma E acting on the solute particles of mass m. Here and in the following greek indices run from 1 to D and we adopt 
Einstein's summation convention over repeated indices. The left-hand side is a conventional streaming operator, while 
the right-hand side is a Fokker-Planck operator with constant friction coefficient 7 and thermal velocity v\ — ksT '/m, 
where ks is the Boltzmann constant and T the temperature of the system. 

In the Lattice-Boltzmann method the continuous velocity v is replaced by a finite set of discrete velocities Vj, 
i = l...b which are vectors on a lattice. Accordingly, the distribution function /(x, v;i) is replaced by b functions 
gi(x; t) tx /(x, Vj; t), and Eq. Q by b equations for each of the gi. In these equations the Vj are no longer variables, 
but fixed parameters. The positions x are discrete points on the lattice, whose size and boundaries are modelled on 
the geometry of the physical problem. Completing the discretization, time t is considered to vary in multiples of a 
discrete step At. This passage from a continuous to a discrete world is schematically pictured in Fig. ^ In order to 
derive the discrete equations, we define more conveniently the external force operator 

C EXT [f] = -a E a d v J (2) 

and rewrite equation as 

(fit + v a d a )f = C[f] (3) 

with C[f] = C FP [f] + C EXT [f\. On the left-hand side we now have a free-particle streaming operator, which can 
be easily discretized, as we will show in the next section. The non-trivial task is to find the correct lattice collision 
operator L corresponding to the continuous operator C. In the BGK case, a systematic procedure has been devised 
in pH Il7| based on Gauss-Hermite quadratures. The procedure is inspired by the pioneering ideas by Grad [lSf 
to solve the Boltzmann equation using the so-called 13-moment system, and has become a useful tool in discrete 
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models of the Boltzmann equation 0, [2(j • It relies on the fact that products of a gaussian with Hermite polynomials 
are eigenfunctions of the BGK operator. We prove in the following that the Gauss-Hermite strategy also allows to 
discretize the Fokker-Planck kinetic equation because BGK and FP operators share the same set of eigenfunctions. 
Indeed they are just different limits of a more general integral operator, devised by Skinner and Wolynes |2l|. The 
additional force term present in C[f] is not diagonal in this basis set, but it has nevertheless a computationally 
convenient form. Although the methodology is not new, its application in the FP context is the novelty of the present 
work. In order to provide a comprehensive and self-contained treatment we give full mathematical details in the 
following. 

The discretization is best carried out by expanding the continuous distribution function /(x, v;t) over a basis set 
of D-dimensional Hermite polynomial tensors 7ia_ (v) (see appendix |A"|> . according to 

/(x,v;^,(v)^ li -4')(x,t)Hf(v) (4) 

where the subscript a is an abbrevation for oti . . . a/, the product denotes contraction on all I indices, and 

e -v 2 /(2vl) 

" (V) = (2^)^ (5) 
is a gaussian weight function, with v 2 = v • v. The expansion coefficients are given by 

F«(x,t)=|dv/(x,v;t)W«(v) (6) 

where the velocity integrals are always taken over M> D . Since 7ia^(v) involve polynomials of order i, the are linear 
combinations of the moments of / 

M^(x,t)= [dvf(x,v;t)v ai ...v am (7) 

with m < I. 

Inserting the expansion in the kinetic equation J2J), we could project the equation on the basis set and derive a 
hierarchy of differential equations for the coefficients. However, that would simply transform the problem into another 
of equivalent complexity. Here we aim instead at a different approach. By means of the Hermite expansion we can 
express the right-hand side of J2J as a function of the moments of /. Using a lattice-discretized distribution function 
we can compute these moments by a suitable quadrature. The expansion involves an infinite number of moments and 
naturally we have to truncate it at a certain order K. We hence assume that 

F®(x,t)=0 iil>K (8) 

and rewrite / as a distribution function that lies entirely in the subspace of Hermite polynomials up to order K 

/(x,v;t) =o;(v)^4rf i l i) ( x '^ ) ( v ) ( 9 ) 

1=0 V T ' 

This assumption is expected to be valid at least for situations close to equilibrium. Using the properties of C FP [f] 
and C EXT [f] detailed in appendix iBl we also find that the outcome of C\f \ = C FP [f] + C EXT [f] lies entirely in this 
subspace. 

The key idea of the Lattice Boltzmann method is that in order to compute the velocity integrals in Eq. (7J , only a 
finite set of velocities is needed. Specifically, we want to compute D-dimensional integrals as discretized sums over a 
fixed set of points. We assume there exists a set of vectors v, € M. D , and a set of real numbers Wi, with i = 1 ... 6, 
such that if p(v) is a polynomial of degree not greater than 2K, the following formula is valid 

b 

/ <2vw(v)p(v) = y^Wip(vj) (10) 



The above equation is then called a quadrature of degree 2K , and the Vj, Wi, the nodes and weights of the quadrature. 
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Because of the truncated expansion Eq. f/u is a polynomial of order K at most. Since we requested a 
quadrature of degree 2K we can now compute the moments of / up to order K, according to 



J dv/(x, v; t)v ai . . . v am = J <2v^j^y/(x, v; t)v ai ...v 



(11) 



b 



\ - /(x,v,;t) \- 

where we defined <7j(x;t) = Wif(jc, v,; t)/u>(vi) and the formula is valid for ?n < .ft*. 

The choice of if is dictated by the application. In practice one is not interested in having / itself, but rather compute 
its moments, which correspond to macroscopic observables. Momentum and energy equations involve moments up 
to second and third order respectively. Consequently it is necessary to require K > 2 or K > 3. The lower order 
moments are labelled by conventional names and the quadratures read 

b 

P = M^= / dv/(x,v;t) = J2 9* ( 12a ) 

J i=l 
b 

J a = pu a = AfW = / dvf(x,v;t)v a = giV ia (12b) 
^ i=l 

P a f 3 = M { ^= dvf(x,v;t)v a vp = ^ giVjgVip (12c) 

i=l 

and if K > 3 we can also compute exactly 

Q a/3l = M^ 7 = / dv/(x,v;i)u o «0U 7 = giViaV^v^ (12d) 

i— 1 

Finding the optimal set {vj, w z } in terms of a minimum number of nodes for a given degree of accuracy, is in general 
an unsolved problem |2^ . However, as far as the solutions of kinetic equations are concerned, it is also important 
that the Vj be vectors of a regular lattice in x space, as shown in Fig. ^ Then for the cases of physical interest 
(D = 1,2,3 and K = 2, 3) a number of possibilities exist with different 6's. The thermal velocity vt can also become 
a free parameter to adjust the quadratures. Some resulting models that are used in practice can be found for example 
in @. 

We now have the prerequisites to create a computational scheme for the Fokker-Planck kinetic equation. Consider 
the distribution function / evaluated at discrete lattice points x and at a finite set of b velocities v, that are also 
lattice vectors. Multiplying both sides of Eq. © by Wi/u>(vi) and using the expansion not on / but on the function 
C[f], we can write 

K 1 

d t g % + v ia d a9i =WiJ2 -WnC^HlHvi) (13) 



where now 

CW = J dvC[f]Hl\v) (14) 
Finite difference time discretization to first-order |17| then leads to 

- 1 

<7i(x + v,At; t + At) - ,g 4 (x; t) = At Wi J2 ~ aTn^a a ( v >) ( 15 ) 

1=0 1 

which defines the lattice Fokker-Planck equation. 

The above expression must be supplemented with an operational expression for the Hermite coefficients Co 1 of 
C[f]. Using the results of appendix iBl they can be expressed as functions of the Hermite coefficient F^p of /. We 
can then write = Ca >,FP + C$' EXT , where 

C W> FP = - 7 /F|) (16a) 
C^ EXT = at FtX +... + alF2-V (16b) 
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are related in turn to the moments Mcp of / by the definition of the Hermite polynomials. 

For K < 3 they read 

F<°> = p (17a) 
= J Q (17b) 
= Pa0-v 2 T pS Q0 (17c) 

^a/3 7 = Q<*Pt ~ v T^af3J 1 + S aJ Jf3 + 8 fa J a ] (17d) 

where p, J a , etc. are related to the gi via the quadratures Eas. l|12|l . 

Putting together Eas. (fT5f) . (fTHfl . and (Pj) . we have then a complete numerical scheme to solve the continuous 
equation J2J). More precisely, we have different schemes according to the choice of the order K in the Hermite 
expansion, which allow corresponding exact computation of the moments of / up to the same order. As an example 
we can write explicitly to second order (K — 2) 

ff^x + ViAtjt + AtJ-ftfx;*) = AtL[g t ] (18a) 

where the lattice collision operator L reads 

E T iV ia Vii3 - Vj. d a/3 

4 



L[ 9l ] = Wl <j [- 7 J Q + + [-2rf(P aP - v^pd af3 ) + a^Jp + ap a ] ia ip 4 1 ap } (18b) 

Urp ZiUri 



Clearly, given the functions gi(x; t) at time t, one can compute the moments p, J ai P a p and hence £[<7i]. The <7i(x; t + 
At) at time t + At are then obtained using the left hand side of Eq. (|18afl . Note that in this way we are not 
calculating the distribution function / but rather g t = Wif(x, v,; i)/w(vi). However, the quantities of interest to be 
sampled are the moments of /, corresponding to hydrodynamic observables. By construction the quadratures provide 
them straightforwardly via Eas. l|12|l . 

The scheme derived here defines an algorithm for the numerical solution of Eq. In the case D = 1, greek 
subscripts are no longer necessary, expressions <|16[) reduce to = —^ylF^ + a^lF' 1-1 ', and Eqs. H17|) simplify as 
well. These expressions coincide with those used in for a second order (K = 2) scheme. In the discrete lattice 
equations are tested against the continuous equation only numerically. In this paper we intend to give a full analysis 
of the scheme defined by Eq. I|18fl ■ Before giving details of the implementation, we check in the following sections how 
reliable the present method is by addressing its stability and the reproducibility of the continuous equation The 
discussion of the computational algorithm is then postponed to sec. 

III. STABILITY ANALYSIS 

In the following we show that Eci. (|15fl can be recast in the linear form 

g'i = J2°»9j (is) 

3 

where q' = q ? ;(x+VjA£; t+At)— ^(x; t) defines a vector g' ', gj = g 3 ;(x, t) a vector g, and dj the so-called collision matrix 
C 0, 123L l24j . For a given lattice geometry, C is a constant matrix which depends only on the operator parameters 
7, cia- In particular we consider isothermal models, where vt is fixed and we make the significant assumption that 
the external acceleration field sl e does not depend self-consistently on the distribution functions gi. The latter case 
will be examined in a subsequent publication. 

The aim of this section is to check for which range of these parameters the scheme embodied in Eq. (|19ll is stable, 
where stability means that upon iterating the scheme, the distribution functions gi (x; t) stay finite at any value of x 
and t. The task is substantially facilitated by the fact that we do not have to first linearize the scheme, as in usual 
Von Neumann stability analysis 8j. By standard arguments in Lattice Boltzmann theory the stability condition 
reads 

|1 + A(<5)|<1 (20) 

where X(C) is any eigenvalue of C. We remark that Eq. 12U|) is a very simple condition valid globally, independently of 
the initial distributions gi (x; 0) or the boundary geometry. Such a feature is an attractive consequence of the linearity 



6 



of the scheme, while much more complicated stability analyses which depend on the local gi(x; 0) are required in the 
full self-consistent LB method |25j . 

Thanks to Eq. H20(l the stability analysis reduces to the spectral analysis of C. We proceed now to first identify 
this matrix, and then compute its spectrum using the results of the previous section. 

As a starting point, we rewrite Eq. (|TK)l as 



n 1 



Nf 
1=1 1 



(21) 



where temporarily in this section we set aside the tensorial notation and enumerate all the terms of the sum (including 
the tensorial contraction) simply from 1 to n. So the index I here represents a shorthand notation for the previous 
set of indices la. Accordingly we redefine the normalization factors as just A 2 , since it is not necessary to know their 
detailed form. We wish then to express Eq. <|21|l in the matrix form of Eq. H19fl using the fact that the dependence 

on gi = Wi /uj(vi)fi is inside C® = J dvd[f\H®(v). 

The first step is then to use the quadratures to write (see appendix lA"|l 



^^wW(v i )W (ro) (v i )=^ OT ^ 2 



(22) 



Defining the matrix Hu = W^(vj) (which contains b rows times n columns), Eq. I|22(l is rewritten in matrix form 

H T WH = A 2 (23) 

where H T is the transpose of H, Wij = WiSij is a b x b diagonal matrix, and A 2 = N 2 Si m is a n x n diagonal matrix. 
Stated otherwise H T (WHN~ 2 ) = I, i.e. WHN~ 2 is a right-inverse {H T )- 1 ' R of H T . 

As a second step, we consider the operator C, and we apply it to / expanded in its Hermite representation 



Cof = Co 



JVi 
i=i 1 



w(v) 



H«(v) 



A 2 



F (i) 



where F"' = / dvfH"\v). Upon projecting along H^ m ^(v) we get 



[ dvH (m) (v)[C'o /] =jr ( ( dvH {m) (v)C< 
^ i=i \J 



w(v) 



A 2 



' i=i 



where the quantity in round brackets defines the elements C m i of a n x n matrix C. 
The third step is to express in terms of the gi 



= f dv/ftWfv) = I dvLu{v)^—H (l) M 
J J w(v) 



^(Vi) 



(v) 

«W(v,) = X;«W (,) (vO 



(24) 



(25) 



(26) 



where the last term can be rewritten in matrix notation as H T g. 

Combining the results obtained in the above three steps we can write 1|21J) in matrix form 



which identifies the collision matrix 



g'/At = WHN~ 2 CH T g 



C/At = WHN~ 2 CH T = (H T )- 1 ' R CH T 



(27) 



(28) 



For clarity we can equivalently write H^ b Cbb = AtC nn H^ b where the matrix dimensions are indicated by explicit 
subscripts. 

Eq. (|28|l is a representation of the collision matrix that allows its spectral analysis. Indeed the spectrum of C is 
directly connected to that of C. In the square case n = b the two matrices are similar and have the same spectrum. 
In the general rectangular case, since b > n, the spectrum of C is contained in the spectrum of C (an eigenvector of 
C being just (H T )~ 1,R v where v is an eigenvector of C). The additional b — n eigenvalues are just 0. 
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We have therefore reduced the problem to the computation of the spectrum of C. We can deduce an explicit 
representation of this matrix using relations ((16(1 and the defining equation ((2511 . The matrix C reads 



/ 



-7 






-7 




-2 7 






-27 



(29) 



where (& E ) contains only of components, the square matrices are diagonal, and all the remaining elements are zero. 
The Fokker-Planck operator fills the diagonal while C EXT occupies the part below. The resulting matrix is triangular 
and the eigenvalues are just given by the diagonal elements, independently of the off-diagonal ones, i.e. 

A fe = k = . . . K (30) 

Consequently, the collision matrix C has eigenvalues A^Ai. Going back to conditions 1)20(1 . the most stringent one 
is for k = K and reads 



< jAt < 2/K (31) 

which in the case of the K = 2 scheme of Eq. I|18|) reduces to < 7 At < 1. These inequalities completely identify 
the range of model parameters for which the scheme proposed in sec. [H] does not lead to an unbounded growth of 
the distribution functions with time. Note that the stability requirement imposes conditions only on the parameter 7 
independently of the external field a. E . This is due to the initial assumption that the field does not depend on the gi. 
Inclusion of self-consistent force fields, that depend for example on the local density Eq. Ill2a[l . would require a more 
careful analysis pfij . 



IV. CHAPMAN-ENSKOG EXPANSION 



A kinetic equation describes a system at the microscopic level of the distribution function /(x, v;t). Define the 
Knudsen number e as the ratio between the mean distance between two successive particle collisions and the char- 
acteristic spatial scale of the system (e.g. radius of an obstacle in a flow). If this number is very small the details 
of particle collisions can be neglected and the system can be considered as a continuum. Using the Knudsen number 
as an expansion parameter, Chapman and Enskog were able to derive from the Boltzmann equation the evolution 
of the hydrodynamic variables (corresponding to the first moments of /) in the continuum limit, thus reproducing 
the macroscopic Navier- Stokes equations pj. Eventually, the expansion has also been used in the context of the LB 
method to derive the macroscopic equations obeyed by Lattice Boltzmann models. The fundamental hydrodynamic 
equations were recovered consistently pfij . 

The Chapman-Enskog procedure is not restricted to the Boltzmann and Lattice-Botzmann equations. In this section 
we apply it to the continuous Fokker-Planck kinetic equation and to the second-order (K — 2) lattice scheme of 
Eq. 1(18(1 . We can then check if the same macroscopic equations for the first moments are reproduced. 

In the continuous case the expansion is straightforward. Indeed one can avoid it completely and obtain the equations 
for the macroscopic variables by just multiplying Eq. Q by v ai . . . v am and integrating over velocity space. In general 
at order m, one obtains the time-derivative of the m-th moment plus the divergence of its flux on the left-hand side. 
On the right hand side the moments of C[f] can be calculated using the Hermite expansion and the properties of 
Hermite polynomials, as was done in the previous section to compute the collision matrix. Explicitly up to order one 
the result is 

d t p + d a J a = (32a) 
d t J a + df}P a0 = -j(J a -pu E ) (32b) 
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where we have introduced the external velocity = a^/j. The first is the continuity equation, the second gives 
the evolution of the first moment J Q , but involves also the unknown second moment P a p. Indeed this procedure 
simply constructs a non-closed hierarchy of equations for the moments of /. However, we are not interested here in 
reproducing the Navier-Stokes equations, nor are we interested in obtaining a closed set of equations. What we wish 
to check in the following is whether the lattice scheme of sec. ITTl actually reproduces the same hierarchy of equations. 

In the discrete case we must make use of the complete Chapman-Enskog expansion. To make it more transparent 
we have divided this derivation into subsections. 



A. Preliminaries 



The macroscopic phenomena that we want to reproduce can occur on different time and spatial scales. For example, 
there may be elastic effects, such as sound propagation, with short time scales, and viscous effects, such as damping, 
with longer time scales. The idea of the Chapman-Enskog expansion is that assuming such a separation of scales, these 
phenomena can be analyzed with multi-scale asymptotic methods |27|. We expand then the populations g^ and the 
spatial and time derivatives in powers of the parameter e, the Knudsen number. The hydrodynamic limit corresponds 
to e < 1. In this limit, noticeable spatial variations take place typically over distances of order er l . Hence, it is 
natural to introduce a macroscopic space variable defined as xi = ex. If we expect to have both propagative and 
diffusive behavior, we must expand up to second order in time, because in diffusion processes inhomogeneities at the 
space scale will relax on the e -2 time scale. Therefore we introduce two time variables t\ = et and t^ = e 2 t. As 
usual in multi-scale methods, we then write 

A = sf+e^+e 2 ^ (33) 
dt = ed\ 1] + e 2 d{ 2) (34) 
d a = ed^ (35) 

Eq. (|33fl defines a corresponding expansion of the moments of g as 

p=5> = E M 0) + 6 * (1) + ^] = p{0) + epil) + eV2) (36) 

i i 

and analogously for J a , P a j3 . For convenience we also rewrite the lattice collision operator Eq. I|18b(l as 



L[9i\ = -jJ a —Wi - 2jP a /3 — -j Wi (37) 



Via _ p. ViaVip - V T d a {3 

where 



4 ' ^ 24 



J a = 3 a - P< (38) 



Pa0 = Pap-V^pSap-^U^Jp+UpJa) (39) 

Since J a and P a p depend linearly on the moments p, J a , P a /3i we can write for them an expansion similarly to Eq. 1)36(1. 
Namely J a = Ja^ + eJa^ + ^Ja^ and analogously for P a p- 



B. Expansion details 



The first step is to apply the expansions defined in the previous subsection to both sides of Eq. 1(18(1 . On the left- 
hand side, we first manipulate the streaming operator as usual in Chapman-Enskog expansions for lattice Boltzmann 
models 8]. Since the scale expansion parameter e is small, the populations vary little from one node the next. We 
can approximate the population <?i(x + v^At; t + At) by its Taylor expansion around <?j(x; t), and write up to second 
order in At: 



5i (x + v,At; t + At) - £/i(x; t) = At 



d t + V ia d a + ~^{dt+ V la d a ){dt + Vipdp) 



(40) 



9 



Using next Eas.l|40|l and (|33135|> . the streaming operator [^,-(x+ Vj At; t + At) — ^(x; t)]/At can be expanded in powers 
of e as: 

order e° : (41) 
order e 1 : [9 t (1) + v ia d^]gl 0) (42) 

order e 2 : fo« + v ta d^ + [d™ + ^(d t (1) + v ia d^)(d^ + (43) 

On the right hand side, in the case of the lattice collision operator the expansion acts order by order on the moments 
and we can write 

L = m+L^+m (44) 

where 

m = - 7 JW V -^ Wi _ 2 7 P W WW -^ S °(> W . (45 ) 

Urp Z/Vrp 

for k = 0, 1 and 2. 

The second step is to equate corresponding orders of the expansion. Thus we obtain to order e° 

= - 7 4 0) ^ - 2 7 P$ ^ -^ t^ (46) 



to order e 1 



Q ^ 24 



(47) 



and to order e 2 the equation can be rewritten more conveniently as 



^W^ 0) + # W, (0) )] = -74 2) S^ - 27^?^=^^ (48) 

The third step is to compute the moment equations associated with Eqs. I|4tj[) - (|47[1 . For the zeroth moment equation 
one can just sum both sides of the equations over i, for the next moments one must first multiply by v 7 , v-yVg, and 
so on. Note that the orders of the velocity moments are not the orders of the e-expansion. For each order in e we 
can compute different moment equations. As we will show shortly, for the purpose of reproducing the macroscopic 
equations l|iS2f> we need up to the second moment equation for orders e° and e 1 and only to the first moment equation 
for order e 2 . The computations are carried using the relations of appendix lAl To order e°, the zeroth moment does 
not give any information, the first and second moment read 

= - 7 J 7 (0) (49a) 

= ( 49b ) 

To order e 1 we find for the zeroth, first and second moments: 

d^pW+d^jW = (50a) 
^4°)+fl«Pi« = - 7 J« (50b) 
bPpW + BPQ®, = -2 7 P« (50c) 
And to order e 2 we only consider the zeroth and first moment equations 

#V 0) + 3 t (1) P (1) + J« - ^> J« = (51a) 
8^ Jf + J« + appg* + ^[9f )(- 7 J( 1 )) + = - 7 J 7 2 ) (51b) 

where we made use of (|50all . (|50bt to derive the first, and of (|50b|) . I|50c|l for the second equation. 
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C. Macroscopic equations 



We can now add up the equations at different orders in e and obtain expanded macroscopic equations for the zeroth 
and first moments of the populations gi. The final step then is to reconstruct the derivative operators from the 
expanded ones. 

For the zeroth moment, we construct e 1 - (|50a|l +e 2 - H51a() (order e° does not add anything in this case) and we get 
e9 t (i y o) + ed^jW + e 2 cfV 0) + e 2 d t (1) p (1 ) + e 2 d^J^ - e 2 ^d« J(D = (52) 
For the first moment, we construct e°- (|49a|l +£ 1 - l|50b|l +e 2 - (|51bj) and we obtain 



e 2 ^[9 t (1) (- 7 J«) + ^)(-2 7 pW)] = -e^JV - eVf (53) 



In these equations both the moments of gi (corresponding to the macroscopic variables) and the differential operators 
are expanded up to order e 2 . We can straightforwardly reconstruct the original quantities using relation i|36|) and 
the analogous ones for the other variables. The spatial derivative is reconstructed in the same spirit noting that 
d a X = ed£\xW + eJfW) = ed^X^ + e 2 d { a 1] X (1 \ where X is any of the moments. In a similar fashion, for 

the time derivatives dtX — ed^X^ + e 2 d^X^ + e 2 d^ X^\ where a term of order e 3 was omitted. Inserting 
Eas. (j49ajl . lj49b|l where necessary, equations l|52|l . (|53|l can then be rewritten as 

d t p + d a J a = ^^<9 Q (J Q - puf) (54a) 

9 t J a + dpPap = ~7(J Q - puf ) 

+ J At dp ^P a p - V^pS a p - i(ltf J/3 + UpJa)^ 

d t ( J a - pu E a ) (54b) 



2 

Interestingly, we find that these equations differ from the continuous equations by one additional term in the first and 
by two terms in the second. All corrections are of order 7 At. 

We can gain more insight in these results by rewriting them in a slightly different way. Let g eq be the solutions of 
L[gi] = 0. From the explicit form (|3*7|l we see that the g\ q satisfy J a = P a p = 0, or equivalently 

J a = pu E a = J e a q (55) 

P af3 = VrpSaf) + -{u®Jp + UpJ a ) = P%*p (56) 

With these definitions we can rewrite Eas. (|54|l as 

d tP + d a J a = l^d a {J a - J?) (57a) 



d t J a + d fj P a p = -i{J a -J e a q ) 

+ 1 Md p (p aP - p; 

7At 



dt(J a -J?) (57b) 

This form shows that for a given value of 7, the closer to equilibrium the system is, the closer the evolution of the 
discrete system is to that of the continuous Fokker-Planck equation. 

Eqs. Q57JI are the final outcome of the Chapman-Enskog expansion of the numerical scheme of Eqs. (|18fl . Together 
with the stability results of sec lllll they complete the analysis of the proposed numerical method. Unfortunately, they 
prove that the scheme does not actually solve the continuous kinetic equation J3J, because of the additional terms 
in Eqs. (|57|l with respect to Eqs. I|32|) . However, as just illustrated, we know explicitly the error made. In the next 
section we exploit this knowledge to build a corrected scheme that is able to solve the continuous equation. 
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V. LATTICE FOKKER-PLANCK ALGORITHM 



Having in mind the results of the previous section we provide here a correct numerical procedure to solve the 
continuous Fokker-Planck kinetic equation Q • 

The results of the Chapman-Enskog expansion suggest that by properly redefining the hydrodynamic variables it 
is possible to recover the correct continuous macroscopic equations. Let 

J* a = (l-^Ja + ^rj (58a) 



Then Eas. (|57|) are rewritten as 



where an effective friction 



2 J 

P*p = (l- 7 At)P a/3 + 7 AtP^ (58b) 



d tP + d a J* = (59a) 

dtJ*+d P* p = -l{J* a -puZ) (59b) 

1 _ 1 At 

7 7 2 



(60) 



is introduced. At the level of the Chapman-Enskog expansion the above equations correspond to the continuous 
Eqs. (|32f> . A similar approach was used in [l2l |2S| by redefinition of velocity in the presence of a forcing term. The 
quantities J* and P*p reduce to J a and P a p in the limit 7Ai — ► 0, Furthermore, if J a = J^ q = pu^ we also have 
J a — J a = Ja Q (and the same for the stress tensor). With walls these equalities do not hold in general. Indeed, the 
boundary conditions set J* = 0, whereas J^ q is non-zero when a field is applied. 

A computational algorithm to solve Q can be divided into two parts, the first initializes the simulation, and the 
second is the dynamical evolution of the gi. 

Initialization. If we perform a simulation using the "bare" definition of the moments in the collision operator, but 
sample the "corrected" moments p, J* and P*@, the latter satisfy the continuous Fokker-Planck equation with 
second-order accuracy, see Eq. I|59fl . but with a rescaled friction 7. Suppose we want to simulate a system with 
a friction 70 . Then in Eq. (|6L)[1 we identify 7 with 70, solve for 7 obtaining 

70 (61) 



70 At 

and use this 7 in the simulation. The external velocity must remain unaffected. So if one wants to apply a 
field one must use in the simulation a field such that 0^/7 = a^/70. The initial conditions are set on 
the starred variables, defined by l|58[) using 7, not 70- 

Simulation loop. Given the set of <?i(x; t) at time t, the gi(x; t + At) at time t + At are found, for each x, by the 
following steps 

1. compute the moments of / using (|11|) . or explicitly (|12fl 

2. compute the Hermite coefficients of /, Ea. (|17l) 

3. compute the Hermite coefficients of C[f], Ea. l|16|) 

4. compute the right-hand side of 115(1 

5. compute the left-hand side of (|15fl 

Then the procedure is repeated. At regular times we can sample the hydrodynamic observables of interest 
corresponding to the moments of /. One has to take care however to sample the starred variables, because those 
are the ones that correctly reproduce the continuous equations. 

An extensive literature is available for the implementation of the Lattice Boltzmann method, where important issues 
such as boundary conditions and large scale code optimization have been investigated in depth Most of the LB 
techniques can be directly extended to the Lattice Fokker-Planck method. For further details we advise then the 
interested reader to more specialized articles, such as the performance studies of [29I l3(ij | . 
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VI. NUMERICAL RESULTS 
A. Numerical limits 

Combining the results of sees. IIIII and IIVI one finds that the second-order scheme of the previous section allows one 
to solve the FP equation for < ^At < 1, independently of the external field a E , and the smaller 7 At, the closer the 
lattice solution will be to the continuous one. Moreover, since in the numerical scheme we use the friction 7 given 
by Eq. l|6*T)l instead of the real 70, the stability condition actually corresponds to < 70 At < 2, so that it seems 
possible to simulate systems with a time step longer than the reciprocal friction. These theoretical findings need some 
numerical back-up. To this purpose we consider here two basic examples in D = 1 for which analytical solutions are 
also available and we compare these solutions with the outcome of simulations in the D1Q3 lattice |6(]. We can then 
check the validity of the proposed scheme and set constraints on the range of parameters. 

In the first example, we consider a system with periodic boundary conditions, initially homogeneous at density po, 
and with no initial velocity. A constant homogeneous field is applied resulting in an external acceleration . From 
the solution of the continuous equations (|32|l . the density does not evolve, while the flux J is uniform and evolves as 

J = (p a S /7o) (1 ~ e-^) (62) 

Direct simulation of the system without the 7 prescription of sec.lVlleads to an exponential solution with a wrong rate. 
Including the prescription and sampling the starred moment J*, one obtains the correct result, as shown in Fig. [21 
However, for 7oAi > 1 (7At > 2/3), we find that the numerical results deviate from the continuous solution, so that, 
even if possible in principle, it is necessary in practice to constrain also the friction 70 to the range < 7oAi < 1. 

In the second example, we consider the same system, but with bounce-back no-slip reflecting boundary conditions 6] . 
A constant field is applied resulting in an external acceleration ajf. Accumulation due to migration results in a 
concentration gradient which is the source of a diffusive flux opposed to the applied field. From the balance of fluxes, 
we find at equilibrium the barometric law for the density 



p(X) (X exp — =-X 

V<4 / 



(63) 



where v\ = hsT/m is the thermal velocity. The same result is obtained from the direct solution of Eas. (|32|l using 
the assumption that the tensor P a p has already relaxed to its equilibrium value P^ given by Eq. (|56|l . Simulations 
without the 7 prescription give an exponential profile, but the exponential slope wrongly shows a dependence on the 
friction 70. With the correct prescription the profile is still exponential, and in order to check the slope, we first 
rewrite the fraction in the right hand side of (|63|l as ciq / (70-D0) where, from Einstein's relation, Dq — v^/jo- From 
an exponential fit of our data we can then derive a simulated diffusion coefficient D dividing 0^/70 by the slope. 
The numerical results are reported in Fig. |21 compared to the continuous value Dq. Slight deviations can be observed, 
especially for small 70. These findings can be explained using the outcome of the Chapman-Enskog analysis. At 
steady state equations become 

d a J* a = (64) 

dpKp = -70£-p«£) (65) 

Using the assumption P*p = P^ = v^pS a p + \ {u^Jp + u^J a ), and and the definition of J*, Eq. Q58[l. we arrive at 
the equations 

d a J* = (66) 
v^p+^—l-^tu^ufdpp + up^] = -^J*-pu%) (67) 

for the redefined variables p, J*. Note that in this case 7 must be identified with 70 above. Using the equilibrium 
result J* = 0, Eq. (|67|) easily yields an exponential solution for p(x) in dimension D = 1, from which we can derive 
the simulated diffusion coefficient D as 

D = ±-** ia »y ( 68) 

7o 2 7 ^ 



The first term is Einstein's relation and the second gives a correction which is small for vanishing external fields. 
The result in Eq. (|68|l is in accordance with the values reported in Fig. |3] since the correction is larger for small 70. 
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Eq. H68fl can also be written as 



7oAt (vF_ 



) 



2" 



D = D 1 




(69) 



where u E — a^/70. Then another way of interpreting the correction is to say that our numerical scheme is more valid 
in a low Mach number regime, i.e. u must be small compared to vt, which numcr ically is 1/a/3 ~ 0.6 for D1Q3 and 
most common lattices. 

Summarizing, we have found that the scheme works, but the theoretical range of parameters must be restrained. 
The physical friction 70 which can be simulated must be such that < 70 At < 1. A small 70 is good to obtain a 
discrete evolution closer to the continuous one, but it must not be too small compared to the external acceleration 
ay since otherwise the low Mach number assumption would fail. As a final remark, in the case of spatially dependent 
forces, this must be true for all the points of the system, as we show in the next section. 



In this section we apply the LFP method to two cases for which no obvious analytical solutions are available. 

In the first case, consider the ID system of the previous subsection with reflecting boundaries and simulations 
on the D1Q3 lattice. We focus here on the time-dependent approach to equilibrium. Such a condition corresponds 
to sedimentation caused by gravity. We also combine the FP collision operator with a BGK operator with a single 
relaxation time r to account for collisions between particles, and possibly hydrodynamic interactions. The stability 
analysis of sec. lllll can be carried straightforwardly also in this combined case, giving the constraint < 27At+At/r < 
2, where 7 At in our scheme is given by Eq. (|61|) . Given a value of 70 in accordance to the previous subsection, we 
can then afford a value of At/r up to 2 — 270 At and slightly above. We report the results in Fig. Interestingly we 
find that the presence of BGK collisions delays the start of the relaxation, but then makes it converge more quickly 
once started. 

In the second case, we consider the 2D system represented in Fig. \5\ where we combine sedimentation and a 
centrifugal force. We apply bounce-back reflecting boundaries on a D2Q9 lattice |(| of 21x41 points. We consider a 
system with friction 70 = 0.1 under the influence of a gravity g = 0.01 and a centrifugal force due to a rotation of 
frequency u r = 0.03. Also at the x borders the low-Mach assumption is satisfied. We report the results in Fig. 
At short times the pure FP system departs earlier from the homogeneous situation. However, at longer times the 
presence of BGK collisions makes the system approach faster the equilibrium distribution. Around t = 400 both 
systems are converged and the final profiles are in agreement with the analytical Boltzmann law. These findings are 
in accordance with the ones of the previous example. 



In order to describe the time evolution of highly asymmetric systems, involving widely different length and time 
scales, like colloidal dispersions, we have extended the Lattice-Boltzmann formalism for the description of fluid flow by 
replacing the standard BGK collision operator by a discretized Fokker-Planck operator to account for the dissipative 
coupling of large solutes to a continuum solvent, without resolving the molecular scale of the latter. Using an 
expansion of the continuous one-particle distribution function in a truncated Gauss-Hcrmite basis, as well as standard 
quadratures with appropriately chosen weights, we were able to reduce the initial continuous Fokker-Planck equation 
to a simple matrix form. The stability of the discrete time evolution is determined by the diagonal elements of the 
triangular collision matrix, which are proportional to the friction coefficient 7. A standard Chapman-Enskog expansion 
leads back to the usual conservation equations derived from the continuous FP equation in the limit At — > 0. For 
finite time steps At, the correct continuum equations are recovered by properly redefining the hydrodynamic variables, 
i.e. by introducing the starred current and stress tensor of Eqs. (|58|l . This leads then to the Lattice Fokker-Planck 
algorithm of sec. 

This algorithm was first tested against known analytical results for the time evolution of simple model systems. 
The numerical efficiency was tested in the non-trivial case of colloid sedimentation in the presence of gravity and 
a centrifugal force. We intend to use the LFP algorithm to investigate ion translocation through heterogeneous 
nanopores, ion transport in swollen clays, and various applications in dissipative colloid dynamics. These applications 
will benefit from the extensive experience gained over the years with the related LB method. 

Due to the well-known mapping between the Fokker-Planck and the imaginary time Schrodinger equation |14| , the 
present LFP scheme is also applicable to the solution of ground-state quantum problems. 



B. Further examples 



VII. CONCLUSION 
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The present LFP scheme has a number of limitations. First of all, the Fokker-Planck equation itself is never fully 
rigorous, since a separation of time scales is never complete, as had already been recognized by H.A. Lorentz [31| so 
that non-mar kovian corrections are always present Secondly, to account for collisions between particles, a BGK 
term may be added to the discrete FP operator, as stressed several times in this paper, and illustrated in two of the 
numerical examples (cfr. Fig. 0] and [BJ). However it is not clear how such a term could account for the long-range 
hydrodynamic interactions between Brownian particles induced by the solvent back-flow. A third limitation emerges 
from the stability analysis, which restricts the range of possible values of the inverse time scales 7 (associated with 
the FP operator) and 1/r (characterizing the BGK operator). Clearly there is a need for an algorithm valid to higher 
order in At. Work to improve along these lines the Lattice-Fokker-Planck method put forward in this paper is in 
progress. 
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APPENDIX A: D-DIMENSIONAL HERMITE POLYNOMIALS 



A complete set of orthonormal polynomials in D variables can be obtained by products of Hermite polynomials in 
a single variable. A detailed presentation can be found in the excellent work of Grad 33|. Here we sketch the basic 
notions and concentrate on the relations which are useful in the present work. 

Consider the space of real functions /(x) of D variables for which the integral J dvw(v)/(v) 2 exists, where w(v) is 

the gaussian weight function defined by Eq. JSJ . A D-dimensional Hermite polynomial of order / is a tensor Tl^ (x) 
of rank I. Each component is a polynomial function in this space. These polynomials form an orthogonal set in the 
sense 



dvw(v)- 



1 



l+m 9L 



(Al) 



where the quantity is zero unless the subscripts a = a\...ai are a permutation of f3 = /3i . . . /3;. It is a sum of 
l\ terms, each one being a product of I Kronecker <5's with the subscripts given by the all possible permutations of 
indices from the two sets a and (3. The first few polynomials read 



W<°>(v 
w£>(v 



1 

VuV/3 



V a VpV 1 - V^(v a Sp 7 + Vp5 al + v 7 5 a p) 
V a VpVjV6 - VT(v a Vf3$~/S + VaV-fSfjg + 

Wy((5 a /3<5 7 5 + SaySps + SasSpj) 



(A2) 
(A3) 
(A4) 
(A5) 
(A6) 



Hermite polynomials form a complete set and the expansion Q is valid, where the coefficients are given by Eq. JBJ. 
In D = 1 the polynomials defined here reduce to the so-called Hermite- Chebyshev polynomials, which differ in 
normalization from the usual Hermite polynomials |34| ■ The derivative of Hermite polynomials satisfies two important 
properties. The first relates an Hermite polynomial of degree / to one of degree I — 1 



The second is the recurrence relation 



pa 



(v)] 



(A7) 



(A8) 
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where (3a denotes the I + 1 indices 0a\ ... a/. 

Making use of the quadratures, Eqs. (jlOJl in sec. [HI integrals of products of Hermite polynomials and a gaussian 
can be rewritten as discrete sums on lattice vectors. The maximum order of the polynomial involved is dictated by 
the order of the quadratures. In the practical case of interest here, a quadrature of order 4 is used in the model of 
Eq. [^J where K = 2. The orthonormality relations Eq. IjAlfl become then the formulae 

^w, = 1 (A9) 

i 

WjVigVip = «4<W (A10) 

i 

22 Wi{v ia Vip - V^Sa^^ViS - V^5 7 s) = Vj^(S al Sj3S + SasSp-t) (All) 
i 

and the remaining combinations, such as J^- WiVi a etc. are simply 0. From the last Eq. (|A11|) we can derive the the 
fourth-order tensor formula 

y^WiVigVipV^ViS = VriSapSyS + S^SpS + SagSpy) (A12) 

i 

APPENDIX B: EIGENFUNCTIONS OF THE D-DIMENSIONAL FOKKER-PLANCK OPERATOR 

The D-dimensional Hermite polynomials defined in the previous section can be used to construct eigenfunctions of 
the Fokker-Planck operator C FP [f] of Eq. £[J in the form of products w(v)h| (v). 
Using the fact that for a gaussian 

d Va u,(v) = -%w(v) (Bl) 

we can write for the action of C FP on these functions 

^ P Kv)«f(v)] = 7^(»fl+^)Mv)«f(v)] 

= yu;(v)(-vpd vp +v*d V0 d V0 )HlHv) (B2) 

Because of relation (|A8|I , we can write 

V$d V0 d Vfs Hg> (v) = DUf (v) + vpd v ,Uf (v) - d Vf ,H^ (v) (B3) 

Using then property (|A7|) it is easy to prove that 

^W^ 1) (v) = -(/ + J D)W«(v) (B4) 

Bringing all the above relations together we get 

C FP [ U (v)«f (v)] = - 7 ^(v)W« (v) (B5) 

which is the eigenvalue property we wanted to prove. From this it is immediate to prove relation (|16afl using (|25|l and 
the orthonormality of the polynomials. 

We can use the above results to prove also relation (|16b(l . From (IBljl and (|A8|) we get 

C BXT Kv)Hf(v)] = -affluXvjHW (v)] 
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from which, using the expansion jlf. 



dvH^{v)C EXT [f] = M« ra, (v)^M^(v)] 



— 1 t?(™-l)„g^( m ) 

= < F ^i + ■■■+ <£KL ( B7 ) 

where we used Eq. <|A1|I and the fact that F ( m ~ 1 ) is invariant under permutations of its m — 1 indices. 
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FIG. 1: The discretization of Lattice Boltzmann methods. A lattice in x space is chosen and the velocity v is restricted to 
a finite set Vj, i = 1 ... b of lattice vectors. The choice of the lattice and the is not arbitrary, but is constructed to allow 
the calculation of the moments of / by quadratures. In the figure a two-dimensional cubic lattice is depicted with 9 velocity 
vectors, usually called the D2Q9 model. 




FIG. 2: A constant field is applied by imposing an external acceleration a® — 0.1. As a consequence a flux J*{i) is created 
in an initially homogeneous system at density po- Because of friction forces, the quantity J* / po saturates at the steady state 
value 0^/70 . Time t is in units of At and the vertical axis has units of Ax/ At where Ax is the lattice spacing. The algorithm 
reproduces the continuous solution (solid lines), but only in the range < 70 At < 1. For 70 At > 1 (inset) discrepancy with 
the analytical solution is found. 
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FIG. 3: A constant external acceleration o is applied in a confined system initially homogeneous. The density converges to the 
barometric law from which a simulated diffusion coefficient D can be extracted. The continuous lines correspond to Einstein's 
relation. The deviation are caused by the discretized solution and can be explained with the results of the Chapman-Enskog 
expansion. The acceleration is in units Ax/ At 2 and the diffusion coefficients in units Ax 2 / At, where Ax is the lattice spacing. 




FIG. 4: Log-linear graphs of the density p(x) as function of the position x at six different times for a ID confined system with 
friction 70 = 0.05 under the influence of an external acceleration ag — 0.01. On a D1Q3 grid of 101 points, the system is 
initially prepared at a homogeneous density po = 1 (horizontal dotted line in the figures) and no initial velocity. As the field 
is applied the system gradually evolves to the barometric law, diagonal dotted line in the figures. We show the evolution for a 
pure FP collision operator (At/r = 0, full line) and combined with a BGK operator (At/r = 1.9, dashed line). 
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FIG. 5: The 2D system. Sedimentation under gravity is combined with a centrifugal force. The external acceleration is a 
combination of the constant g and the linear escape term oj^x. 
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FIG. 6: Contour plots of the density p(x, y) at four different times. The system is initially homogeneous at density po = 1 
and no initial velocity. As the field is applied the density gradually evolves to the characteristic parabolic profiles of the 
equilibrium distribution. We show the evolution for a pure FP collision operator (Ai/r = 0) and combined with a BGK 
operator (Ai/r = 1.8). At t = 400 both systems are converged to the analytical Boltzmann law. 



